function [V_coord] = ic_gen(fname,nx,ny,lx,ly,nv)
%function ic_gen(fname,nx,ny,lx,ly,nv) generates a collection of nv
%periodically copied Taylor vortices on the grid with (nx,ny) gridpoints in
%a box of size -lx/2 < x < lx/2 and -ly/2 < y < ly/2. The parameters of
%vortices that are generated are random.
%
%The result is written to output file fname in format readable by C code. 

global X Y
dx = lx/nx;
dy = ly/ny;
x = -0.5*lx + (0.5:nx-0.5)*dx;
y = -0.5*ly + (0.5:ny-0.5)*dy;

[X Y]=meshgrid(x,y);

w = 0;
V_coord = [];
for i = 1:nv
    
%    x0 = 0; y0 = 0; U = 1; a = 0.25;
    x0 = lx*(rand - 0.5);
    y0 = ly*(rand - 0.5);
     
    U = 0.8*rand + 0.2;
    a = 0.7*rand + 0.3;

    w = w + taylor_vortex(a,U, x0,    y0) + taylor_vortex(a,U, x0-lx, y0   ) + ...
            taylor_vortex(a,U, x0+lx, y0) + taylor_vortex(a,U, x0   , y0-ly) + ...
            taylor_vortex(a,U, x0,    y0+ly) + taylor_vortex(a,U, x0-lx, y0-ly) + ...
            taylor_vortex(a,U, x0+lx, y0-ly) + taylor_vortex(a,U, x0-lx, y0+ly) + ...
            taylor_vortex(a,U, x0+lx, y0+ly); 
    V_coord = [V_coord; a U x0 y0];    
end
w = 2*sin(X+Y);
figure(1)
surf(X,Y,w)
shading interp
w=reshape(w,nx*ny,1);
file_id=fopen(fname,'w');
fprintf(file_id,'%e\n',w);
end

function out = taylor_vortex(a,U,x0,y0)
  global X Y
  dX = X - x0; 
  dY = Y - y0;
  ov2 = 1./(a*a);
  out = (U/a)*(2-(dX.^2+dY.^2)*ov2).*exp((1-(dX.^2+dY.^2)*ov2)/2);

end


